Sensors 2014, 14, 9429-9450; doi:10.3390/sl40609429 



OPEN ACCESS 



sensors 

ISSN 1424-8220 

www.mdpi.com/journal/sensors 

Article 

Arterial Mechanical Motion Estimation Based on a Semi-Rigid 
Body Deformation Approach 

Pablo Guzman '*, Ghassan Hamarneh 2 , Rafael Ros 3 and Eduardo Ros 1 

1 Department of Computer Architecture and Technology, ETSI Informatica y de Telecomunicacion, 
CITIC-UGR, University of Granada, 18071 Granada, Spain; E-Mail: eros@ugr.es 

2 School of Computing Science, Simon Fraser University, Burnaby, BC V5A 1S6, Canada; 
E-Mail: hamarneh® sfu.ca 

3 Hospital Universitario San Cecilio, Servicio de Angiologfa y Cirugfa Vascular, 1 8008 Granada, 
Spain; E-Mail: rafael.ros.sspa@juntadeandalucia.es 

* Author to whom correspondence should be addressed; E-Mail: pguzman@ugr.es; 
Tel.: +34-958-24-6128; Fax: +34-958-24-8993. 

Received: 31 January 2014; in revised form: 18 April 2014 / Accepted: 21 May 2014/ 
Published: 27 May 2014 



Abstract: Arterial motion estimation in ultrasound (US) sequences is a hard task due to 
noise and discontinuities in the signal derived from US artifacts. Characterizing the 
mechanical properties of the artery is a promising novel imaging technique to diagnose 
various cardiovascular pathologies and a new way of obtaining relevant clinical information, 
such as determining the absence of dicrotic peak, estimating the Augmentation Index (AIx), 
the arterial pressure or the arterial stiffness. One of the advantages of using US imaging is 
the non-invasive nature of the technique unlike Intra Vascular Ultra Sound (IVUS) or 
angiography invasive techniques, plus the relative low cost of the US units. In this paper, 
we propose a semi rigid deformable method based on Soft Bodies dynamics realized by a 
hybrid motion approach based on cross-correlation and optical flow methods to quantify 
the elasticity of the artery. We evaluate and compare different techniques (for instance 
optical flow methods) on which our approach is based. The goal of this comparative study 
is to identify the best model to be used and the impact of the accuracy of these different 
stages in the proposed method. To this end, an exhaustive assessment has been conducted 
in order to decide which model is the most appropriate for registering the variation of the 
arterial diameter over time. Our experiments involved a total of 1620 evaluations within 
nine simulated sequences of 84 frames each and the estimation of four error metrics. We 
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conclude that our proposed approach obtains approximately 2.5 times higher accuracy than 
conventional state-of-the-art techniques. 

Keywords: computer vision; ultrasound; wall motion; arterial stiffness; elastography; 
carotid; motion analysis 



1. Introduction 

Estimating the variation of motion in the artery for vascular characterization [1] is a new technique 
that helps doctors to detect specific diseases. Other non-invasive techniques such as Ankle Brachial 
Pressure Index (ABPI) [2] or Augmentation Index (AIx) have been used to estimate parameters (blood 
pressure) that are associated with peripheral vascular diseases. For example, Mortensen et al. [3] 
demonstrated the relation of AIx and the Marfan syndrome, the role that involves AIx in the 
hypertension field [4] and the increase of the arterial stiffness in human subjects with Type 1 diabetes 
mellitus [5]. The way to estimate the pressure parameters becomes very limited due to the fact that 
such measures cannot be estimated in other parts of the body besides the carotid artery where we have 
an easy access with US. Arterial pressure and arterial wall motion are related since estimating the 
pressure requires measuring the variation of the diameter of the artery, as it is indicated in 
Equation (2). The importance of the wall motion artery's characterization has been also discussed by 
several authors who have demonstrated that radial [1,6-8] and longitudinal [9-11] motion are 
promising indicators to be associated with certain diseases or pathologies. 

Existing commercial solutions such as Tissue Doppler Imaging (TDI) focus on the velocity 
measurement of the myocardial motion using Doppler principles. This technique has been extended to 
other applications in echocardiography [1,6], to determine the mechanical properties of vessels by 
means of TDI. The main problems with using TDI are that the motion vector measurement can only be 
done in parallel to the direction of the ultrasound beam, TDI measures absolute tissue velocity, and it 
is not able to distinguish all passive motion [12]. 

Different solutions have been proposed to characterize the wall artery motion directly from 
ultrasound images in order to complement the information about motion patterns extracted from 
B-mode US. Image intensity correlation techniques have been widely used in ultrasound due to their 
robustness under noisy environments. Golemati et al. [13] compared the displacement error produced 
in block matching [14] and optical flow [15,16] methods over a simulated dataset. The matching 
feature is also an important factor, where Soleimani [17] demonstrates that by including the gradient in 
the local search, the method results improve. The inclusion of a Kalman filter [18] to update the 
reference block and the displacement vector [19-21] has been also evaluated. This method becomes 
useful when the registered data is corrupted by significant amounts of noise, but in cases where the 
information is not corrupted at all, the filter does not improve the accuracy of the system rather 
produces over-smoothing. Other authors [22,23] go one step further and not only measure the 
displacement of the wall, but also include the Pulse Wave Velocity (PWV) to estimate the pressure by 
mean Moens-Korteweg Equation (1), that relates the PWV to the elasticity of the arterial wall: 
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PWV = 



p dA 



(1) 



where A Q is the arterial diameter in the diastole, p is the density of blood, dA is the variation of the 
diameter over time (determined using cross-correlation in [22,23]) with respect to the artery in resting 
position, and dP is the difference in the pressure with respect to the end-diastolic pressure P(0). The 
pressure produced in the artery can be estimated using PWV, as shown in Equation (2): 



Compared to previous works, our method supports sub-pixel accuracy and incorporate collective 
motion information to define the wall artery motion. The major contribution of this paper is the 
evaluation of different methods and how they can be integrated to better address our problem of 
estimating the change in diameter dA. In this work, in order to enhance existing motion tracking 
methods, a combination of similarity transformation, non-rigid deformations, statistical filtering, and 
hybrid motion estimation techniques are proposed. In this way, it will be possible to estimate useful 
parameters instead of using more expensive and invasive methods that put the patient's well-being 
at risk. 

The paper is organized as follows: Section 2 introduces a detailed explanation of the evaluated 
models, the process to generate the ground truth estimation and the combination of the methods that 
will be evaluated and compared in Section 3. In Section 4, the obtained results of the evaluated 
methods will be discussed and finally, Section 5 summarizes some conclusions and outlines for 
proposed future work. 

2. Material and Methods 

In this section, the methodological 'building blocks' used in this paper will be first briefly 
described. Then, the performance of the methods in different analysis pipelines will be evaluated. 

2.1. Evaluated Methods 

2.1.1. Block Matching 

The block matching (BM) technique has been a very popular method in the ultrasound field because 
it provides a robust estimation of the motion by means of comparing the similarity between blocks of 
different images. One of the uses of motion estimation via BM technique is the one proposed by 
Basarab [24], where the elastography map is estimated to show hidden objects such as cysts or cancer 
tumors in ultrasound imaging. This work uses a multiscale scheme to avoid errors in the motion 
estimation and to obtain a low sub-pixel resolution. It is important to remark that normalized cross 
correlation (NCC) block matching method is one of the most popular techniques utilized in ultrasound 
tracking [9-11,22,25]. In our evaluation, it was decided to make use of Lewis method [26] due to the 
fact that the obtained performance is much superior to the original one (approximately 15 times faster). 
Lewis method consists of a modification of the NCC technique where the similarity is given by 
Equation (3): 



dP = 



p ■ PWV 2 



dA 



(2) 
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Y.{i,j)ew(h(\.j) - h) ■ Q 2 (x + t,y +;) - h) 



(3) 



where /i is the reference block with size W, h is the image where the correlation is carried out, 7 X and 
I 2 are their respective means. The numerator in Equation (3) can be sees as a convolution (as in 
Equation (4)) between two images, f(i,j) = — h an d t(i,j) = I 2 (i,j) — I 2 , and can be 

efficiently solved by means of a Fourier convolution: 



I 

(i,j)eW 



f(i,j) ■ t(x + i,y +j) 



(4) 



On the other hand, the denominator of Equation (3) must be solved efficiently. Lewis [26] proposed 
the use of the integral of the image technique to compute it efficiently and reduce the computation cost. 

2.1.2. Optical Flow 

The temporal variation in an ordered sequence of images allows the estimation of the optical flow 
2D vector, usually denoted as v = (u, v) , and is computed based on the constant-brightness 
hypothesis, which assumes that the pixel brightness remains constant over time. This leads to the 
formulation of the famous optical flow constraint Equation (5): 

uf x +vf y + f t = 0 (5) 

where u and v are the optical flow components and the spatio-temporal derivates are represented by f x , 
f y and fi respectively. It is important to remark that in the two considered optical flow methods 
(described in Sections 2.1.3 and 2.1.4), the texture domain of the image was used, as proposed by 
Wedel et al. [27], so as to avoid problems with the brightness assumption. 

2.1.3. Optical Flow: Lucas and Kanade 

On the basis of the optical flow constraint equation, Lucas and Kanade [15] proposed the 
minimization of the error Equation (6) using the sum of the least squares: 



E(u, v) = ^V*(i)u + f y (i)v + f t (0y 



(6) 



165 



Equation (6) is minimized by taking the partial derivatives with respect to the optical flow vector v. 
The resulting v minimizes the differential error between the previous image and the current image and 
is given is presented in Equation (7): 

[-£/x(0/t(0" 



V = 



^V(0 £/x(0/y(0' 
iEB 165 

^V*(0/y(0 ^Vy 2 (0 



166 



iEB 



iEB 



iEB 



(7) 



where v is the optical flow vector, the spatio-temporal derivates are represented by f x , f y and f t 
respectively, and the subscript i is the i-th element of the integration block B. 
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2.1.4. Optical Flow: Anisotropic TV-LI 

The anisotropic optical flow Equation (8), proposed by Werlberger et al. [28], is an extension of the 
popular method TV-LI optical flow [29], which is based on a regularized propagation technique 
similar to the one proposed by Horn and Schunck [16]. Werlberger proposes an anisotropic diffusor 
that does not propagate values through the edges, with better preserves image structure. In this case, 
the original formulation was changed in order to reduce the computation cost in the last term of 
Equation (8), where possible artifacts (e.g., occlusions) can be rectified over the time incorporating 
feedback from previous optical flow: 



min U j sup Wd]sl jj } \ [D^u d ) ■ p d - e -|— + — (u d - v d ) 2 




' + 8\p(v(x))\ 



(8) 



where D li = exp (— a|V/|^)nn 7 ' + n^n^ 7 ), n = ^ as the normal vector, and n 1 " the tangent vector 
of a given point. u d is the optical flow vector and u' d is the previous warped optical flow [30,31]. 

2.1.5. Kalman Filter 

In noisy systems, the Kalman Filter [18] has been proposed due to its robustness and efficiency. 
This method is based on a statistical approach to determine the current estimation of a linear system 
from a collection of previous observations over time as described in Equation (9): 

x k = Ax k ^ + Bu k _ x + w fc _! (9) 

where x k is the predicted estimation in the current time k, x k -^ is the previous observation, A and B 
describe the transition and control matrix respectably, u fc _ x is the control signal, and w k ^ t , the process 
noise of the system. 

2.1.6. Similarity Transformation 

Incorporating shape prior knowledge has become common practice in segmentation methods in the 
last decades. Cootes et al. [32] proposes a statistical method able to deform a contour by means of 
weighting relevant eigenvectors (P) by shape parameters (b) with the objective of adapting the contour 
to a desired object in the image, as shown in Equation (10): 

x = x + P-b (10) 

The objective of this method is estimating the shape parameters as well as the pose parameters 
(translation in x and y-axis, scale, rotation) that locate the desired object to be segmented in 
Equation (11): 

v = {t x ,t y ,s,B,b} (11) 

These parameters are found by means of matching each landmark with the previous trained 
normalized gradient profiles and solving linear equations [33]. In our work, we will not adopt the 
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statistically-based deformations, thus the local deformation terms can be discarded. In other words, we 
neglect the shape variability encoded in b and focus on estimating the remaining pose parameters only. 
Therefore, Equation (10) is no longer needed in our approach and we only utilize the weighted 
similarity transformation proposed in [33]. Such a similarity transformation is obtained by means of 
the weighted sum minimization in Equation (12): 

E = (xl-M{s,6)[x2\-t) T W(xl-M{s,6)[x2\-t) (12) 

where: 



ks, v) [ y \ y, . sinB , . x + , g . CQsB , . y 



(s ■ sinO) ■ x + (s ■ cosO) ■ yj 

t (j-x> ty> • • • > tx> ^y) 

xl is the origin point and x2 the translated point, s the scale, 0 the rotation, t the translation vector, 
and W is a diagonal matrix of weights for each point. 

2.1.7. Soft Body Dynamics 

Soft body models have been widely used in computer science to carry out realistic physical 
simulations of motion and deformable objects. Rather than a statistically-based deformation model in 
Equation (10), this paper will focus on the popular mass-spring model, which is based on a mesh of 
nodes (masses) and connected by means of elastic links (springs). The basis of this method relies on 
Hooke's law, to simulate the spring force, and the second Newton's law to simulate the dynamics by 
time integration. In this work, a simplification of the idea proposed by Hamarneh et al. [33] will be 
adopted. The authors describe a system in Equation (14) that involves forces generated by the springs' 
system, in Equation (15), in a controlled environment (Equation (16)) allowing for speedup/slowdown 
of the velocity of the mesh's nodes: 

j _ jHooke _|_ jViscous (14) 

fT~- -*.(lh-^ll-n/)i7^-(^(*.-*/i^)ii7 Z f» (is 

Ay || \ || Aj -*7||/ || i j \\ 

^Viscous = _ kvV , (16) 

where f t is the final estimated force, k s the Hooke's spring constant, x £ the Cartesian coordinate of the 
i-th node, r t j the rest length associated to a link between two nodes, k d the damping constant, k v the 
viscosity coefficient, and v t as the velocity of the i-th node. Once the nodes' force is obtained, it is 
possible to estimate the acceleration, velocity, and position of each node by means of the iterative 
scheme in Equation (17): 

A 

771; 

_ ..old , „ At (17) 



Vi = v ° La + a .At 
%i = x f ld + ViAt 
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where a t is the acceleration, f t the force described in Equation (14), v t the velocity, x £ the position, 
and At the time interval. 

2.2. Ground Truth Estimation via Ultrasound Simulation 

To evaluate the proposed methods, the Field II Ultrasound MATLAB library [34] was used to 
generate nine sequences simulating the wall displacement of the common carotid artery. Each 
sequence involves a complete cycle of the cardiac system with a frequency of 25 Hz per cycle. These 
simulations were generated with 1024 physical elements, a transducer center frequency of 
5 MHz, 100 MHz of sampling frequency, and 64 active elements. The sequences are based on three 
different topologies as shown in Figure 1, where diverse amplitudes of motion are used, as explained 
later in this section. Although these synthetic images are clearer than real US artery images, having the 
ground-truth motion allows quantitative evaluation metrics to compare different methods. We include 
also the simulation parameters used to produce these synthetic sequences (Table 1) to facilitate the 
reproduction of our results. 

Figure 1. Ultrasound images used in the data set test, generated by Field II U.S. simulator. 




Following the steps of Stoitsis et al. [35], that describe a mathematical mechanical deformation 
model of the arterial wall as a separable method in space and time, this paper will simulate the radial 
and longitudinal displacement of the artery by means of Equations (18) and (19), which is a 
simplification of the original method. 

rt(t) = Y\(t 0 , tx) ■ sin 2 ^ + Y\(fi, t 2 )-(a + b- t) (18) 
Y\(ti, tj) = f ■ (1 + tanh(d ■ (t - tj))) ■ (1 + tanh(d ■ (t, - t))) (19) 

where a and / determinate the amplitude of the waveform, b defines the slope in the second part of the 
curve, c and T are coefficients that control the initial part of the curve, d determines the wall artery 
speed, t x and t 2 correspond to the duration of the first and second pulse of the waveform and t is the 
time variable. In the generated data set, the chosen parameters are shown in Table 1, where different 
values of / let us control the amplitude of the artery displacement. For each generated artery topology, 
three values of / will be applied to produce different radial and longitudinal displacements. 
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Table 1. Wall Displacement Simulation Parameters. 

Parameter Value 

a 15.14 

b -0.64 

c 1.5 

t x 0.25t 

t 2 0.65t 

d 1.22 

/ 0.06,0.12,0.25 

r l 



2.3. Proposed Methods 

Following the introduction of the methodological building blocks in the previous section, the 
proposed combination of methods to be evaluated in Section 3 will be explained in detail. To estimate 
the motion of the wall of the artery, two optical flow methods and the block matching approach 
proposed by Lewis [26] will be evaluated. 

One of the main problems in classical block matching techniques is the sub-pixel accuracy. To 
handle this problem in an efficient way, other authors [36,37] proposed a combination of optical flow 
and block matching to increase the motion vector precision. This approach relies on estimating firstly 
the motion vector by means of the block matching technique, and then applying a warping [30,31] to 
the block and computing the optical flow to estimate the sub-pixel information (as shown in Figure 2). 

Figure 2. Block matching with sub-pixel accuracy by means of the optical flow scheme. 
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Most of the wall artery tracking papers take into account the tracking individually [13,20,21] but do 
not consider all the tracking points as a set of data that define a semi-rigid object in motion. In this 
paper we use all tracked points when estimating the similarity transformation parameters (rotation, 
scale, and translation), as explained in Section 2.1. This method obtains all the motion vectors for all 
point and estimates the transformation parameters by means of Equation (13). One of the advantages 
of this method is the robustness against noise, allowing the computation of the parameters that define 
the transformation even with some wrong (or outlier) motion vectors. Figure 3 describes the steps that 
define this proposed scheme. 

Figure 3. Similarity transformation given the motion vectors obtained with the hybrid 
BM-optical flow method. 
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When the similarity transformation is applied, it produces an error in the position transformation. 
This artifact happened because the method considers the systolic-diastolic-systolic transition as a 
scale-translation transformation over time. To reduce this error, rather than using the statistical shape 
deformation model of Equation (10), a physically-based simulation method (Soft Body, Section 2.1) 
could be included to the upper and lower set of points (individually) that track the artery, as indicated 
in Figure 4a. The pipeline of this approach is illustrated in Figure 4b. 

Figure 4. (a) Pipeline of the proposed method with physics simulation (mass-spring) and 
(b) Illustration of the spring connections in an ultrasound image. 
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Finally, the KF (Kalman Filter, Section 2.1) will be incorporated into this evaluation. This filter was 
previously used by Gastounioti et al. [20,21] to increase the accuracy of tracking the wall artery 
tracking. In this particular case, the KF will be incorporated into different proposed schemes as 
described by Figure 5 with the objective of being evaluated a posteriori in Section 3 as well as the 
other proposed schemes. It is important to note that in our evaluations a simple updating scheme was 
added to update the reference block (the initial tracked block) in each frame of the sequence, as shown 
in Equation (20): 

B-newref = a ' B 0 ld_ref + (1 _ <*) " B estimated (20) 

where B new _ re f is the new estimated reference block, B old _ re f the old referenced block, B estimated is 
the estimated displaced block, and a is the parameter that controls the amount of information that 
remains from the old reference block; we empirically set a to 0.98. 
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Figure 5. Incorporation of Kalman Filter in the previously proposed schemes. Each 
scheme is encoded with a different color. 
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3. Results 



As described in Section 2.2, we generated a set of nine sequences that simulate the wall artery 
motion with different amplitudes and topologies. The displacement of the wall of the artery, in our 
evaluation set, varies from 2.50 up to 9 pixels. To evaluate the proposed methods in Section 2.3, a 
similar metric as the one used by Golemati et al. [13] will be used. In our evaluation, each sequence 
contained 84 frames and involved three cardiac cycles. To obtain an impartial evaluation, six different 
positions (a total of 54 evaluations per method) will be evaluated in each sequence and the Cartesian 
coordinate error Equations (21)-(23) and the diameter error over the time Equation (24) will be 
measured by mean of root-mean- squared error (RMSE): 



M-1 W-1 



-long 



N -M 



j=0 i=0 



M-1 W-1 



^rad 



In 



(21) 



(22) 



J = 0 1 = 0 




M-1 W-1 



^ZZ (||pa; ' )_p ' a; ' )||2) 



(23) 



J=0 1=0 



w 

M-1 2 



-1 



N 



— V V(d(ij)-d'(i,;)) 2 



(24) 



where x(ij'), y(i,j) and p(i,j) are the i-th x, y and (x,y) coordinates in the j-th frame of the sequence 
with its respective ground truth x'(i,y), y'(ij'), p'iUj) and is the estimated diameter of the 

artery and its ground truth d'(i,j). N is the number of points, in our case, six points, and M is the 
number of frames per sequence. Our evaluation uses five motion models: Lucas and Kanade, 
Anisotropic Huber LI, BM (Lewis method), BM + Lucas & Kanade, and BM + Anisotropic Huber 
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TV-LI. For each motion model, the methods described in Section 2.3 will be used. These methods 
include the similarity transformation (ST), the Kalman filter (KF), and the mass-spring (MS) physics 
based model. The motion methods will also be evaluated individually as shown in Table 2. 

Table 2. Evaluated models where the used methods are indicated. Not all the combinations 
have been used because some of them were nonsensical. 
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Figure 6. Longitudinal error results after being evaluated with different methods. 
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At first, the estimated longitudinal error will be evaluated by means of Equation (21), in our 
evaluation set (Section 2.2). As can be appreciated in Figure 6, the methods that make use of the Lucas 
and Kanade algorithm generate the highest number of errors. This is due to the fact that this method 
does not handle properly the aperture problem compared to the other motion evaluated methods. Ml 
with block matching and M5 method with anisotropic TV-LI motion estimation produce the best 
results in the longitudinal motion estimation. In some plots the errors values are above the maximum 
values of the plot. We have reduced the plot range to better discriminate among the other approaches. 

In the radial motion evaluation (Figure 7) computed by means of Equation (22), the methods with 
Lucas and Kanade obtain again worse results. But unlike the longitudinal evaluation, the mix of block 
matching and anisotropic TV-LI obtain the best results and a lower deviation with respect to the best 
results of the block matching approaches, achieving almost 50% less error. It is important to remark 
that in both evaluations (longitudinal and radial), the inclusion of the mass-spring method 
(M5-M6) helps, in general, to reduce the position error. 
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Figure 7. Radial error results after being evaluated with different methods. 
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Figure 8. Position error results after being evaluated with different methods. 
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A global position evaluation (Figure 8) is carried out by Equation (23), where the radial and 
longitudinal displacement is taken into account. In general, the block matching (Ml version) and 
anisotropic TV-LI methods achieve the best results with the difference that BM generates less 
standard deviation in the error among different evaluated sequences. 

It is important to obtain a method capable of achieving good position precision. In our case, it is not 
critical to estimate the evolution of the diameter of the artery over time with a model that generates a 
small deviation on the elastomer's position (+1 pixel). To evaluate the most significant parameter, the 
diameter of the elastomer over time, estimations obtained from Equation (24) will be evaluated in our 
assessment set (Figure 9). The results reveal that the best method is the combination of block matching 
and anisotropic TV-LI (M5). If we compare it with the previous results that evaluate the position error, 
the best solution (BM-M1) generates 3.1 times higher error than the new best solution obtained. 
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The combination of optical flow and block matching obtain almost two times higher precision than the 
methods working individually. 

Figure 9. Diameter error results after being evaluated with different methods. 
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After evaluating the different proposed approaches, it can be observed that the inclusion of the 
Kalman filter (M2, M4 and M6) provides worse results with respect to the other approaches. The 
methods that include the similarity transformation reduce the error up to 14% in some proposed 
approaches. After incorporating the soft body model, an increase of 13% in the precision is obtained in 
the last method (BM and anisotropic TV-LI). Finally, it can be concluded that a combination of optical 
flow and block matching and the M5 scheme becomes the most precise technique to estimate the 
desired parameter. 

To show a more detailed evaluation between the best method and other results, a Bland- Altman plot 
is produced. Figure 10 shows the difference between the best method with other evaluated methods 
and their averages. The middle line indicates the average difference of both methods, whereas the 
upper and lower lines represent 95% limits of agreement with 15.80% (Figure 10a), 49.78% 
(Figure 10b), 33.66% (Figure 10c), and 78.39% (Figure lOd) of window (defined by 1.96 times the 
standard deviation with respect to the mean difference) displacement with respect to the origin 
coordinate. It can be concluded from Figure 10 that there is an overall good agreement of the 
amplitudes between the BM and anisotropic TV-LI (M5) method and the reference ones. To facilitate 
the evaluation, further tabulated results are listed in the appendix. In the next section, the obtained 
results will be discussed and the method will be evaluated in real cases with different human subjects 
in vivo, with the objective of validating this technique. 
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Figure 10. Bland-Altman figure, where the best obtained results block matching & 
anisotropic TV-LI (M5) are compared with (a) block matching (Ml), (b) block matching 
(M5), (c) anisotropic TV-LI (M3), and (d) anisotropic TV-LI (M5). 
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4. Discussion 



In this section, the proposed approaches and the impact of each method in different evaluations will 
be discussed. Figures 11 and 12 show the proposed evaluations in two different sequences, in order to 
illustrate the response of the methods under different motion variations. 

Figure 11 corresponds to a simulation that produces a maximum displacement of 18 pixels (the 
variation in displacement between the upper and lower arterial wall), while Figure 12 generates a 
maximum wall displacement of 5 pixels. As was discussed in the previous section, the Lucas and 
Kanade method is not the most appropriate algorithm to register the motion in US images as illustrated 
in Figures lla,b and 12a,b. When displacement vectors are long, the optical flow based techniques are 
not the most convenient ones, because these methods have a maximum limit to determine the motion 



Sensors 2014, 14 



9443 



vector (thus improving the working range would require multiscale schemes such as [38]) as illustrated 
in Figure lid. 

Figure 11. Diameter evolution over time evaluated in a sequence with long displacements 
using (a) Lucas-Kanade, (b) block matching & Lucas-Kanade, (c) block matching, 
(d) anisotropic TV-LI, and (e) block matching & anisotropic TV-LI. 
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Figure 12. Diameter evolution over time evaluated in a sequence with small displacements 
evaluated with (a) Lucas-Kanade, (b) Block Matching & Lucas-Kanade, (c) Block 
Matching, (d) Anisotropic TV-LI, and (e) Block Matching & Anisotropic TV-LI. 
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Figure 13. Diameter evolution of the common carotid artery (CCA) in real ultrasound data 
in different subjects where first row (a) correspond to a healthy patient and the last two 
rows (b-c) belong to patients with presence of atheroma plaques. 
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Moreover, Block Matching techniques do not have this limitation, but produce rough results and 
allow no sub-pixel precision (Figures 11 and 12c). A good solution is the combination of correlation 
and optical flow techniques to avoid this displacement limitation and obtain sub-pixel precision 
(Figures 11 and 12e), acquiring twice more precision than with only block matching method (approx. 
0.25 pixel error), but with the inconvenience of increased computation time. The inclusion of Kalman 
filters does not significantly increase or decrease the results, but it is interesting to include it in 
hypothetical cases when the system has a severe disturbance and noise. The main problem with the use 
of this filter is that, depending on the settings of the parameters, the signal may be over smoothed and 
shifted in relation with the desired one. 

At this point, the proposed methods have been evaluated on synthetic US B-Mode imaging. To 
demonstrate and validate that the best approach is able to work in real sequences, it is evaluated in 
different subjects in the common carotid artery as shown in Figure 13. In this brief evaluation, a 
healthy patient (Figure 13a) was involved, whose diameter motion curves showed clearly the dicrotic 
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peak (attributed to the elastic recoil of the aorta and aortic valve) while the other two patients (with 
presence of atheroma plaque) the dicrotic peak is absent (Figure 13b,c). 

5. Conclusions 

The objective of this work was the design, evaluation and comparison of methods able to 
characterize arterial wall motion. A set of methods has been evaluated with the objective of 
determining which approach better handles our problem, the estimation of the diameter of the artery 
over time. The motion methods were selected according to the obtained results in other works, with the 
goal of comparing our approach against these other methods and evaluating its accuracy. It has been 
demonstrated that our proposed combination of methods based on similarity transformation, non-rigid 
deformations, statistical filtering, and hybrid motion estimation techniques enhance existing state of 
the art approaches, up to 2.5 times more accurate than state of art techniques. 

Synthetic US sequences with different patterns of motion were generated to allow quantitative 
comparative analysis of different methods and combination of techniques. Our experiments involve a 
total of 1620 evaluations within nine simulated sequences of 84 frames each and four error metrics. In 
fact, the assessment that appropriate integration of different techniques has a clear impact on the 
final performance represents an important contribution of this work. Another advantage that must be 
remarked is that the proposed methods supports large displacement vectors unlike optical flow 
techniques that are limited in working range and require multiscale schemes. 
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Appendix: Results 

Table Al. Longitudinal error results after being evaluated with different methods. 

Methods Ml M2 M3 M4 M5 M6 

Lucas &Kanade 62.87 16.95 59.98 73.91 75.99 42.52 

BM+ Lucas & Kanade 52.79 45.05 57.41 39.51 29.76 19.48 

BM 1.10 1.15 1.68 1.68 1.27 1.28 

Anisotropic Huber-Ll 1.43 1.42 1.56 1.56 1.12 1.12 

BM + Anisotropic Huber-Ll 1.73 1.62 1.88 1.87 1.52 1.51 



Table A2. Radial error results after being evaluated with different methods. 



Methods 


Ml 


M2 


M3 


M4 


M5 


M6 


Lucas & Kanade 


7.89 


2.37 


4.96 


8.50 


4.93 


3.33 


BM+ Lucas & Kanade 


4.69 


4.42 


3.96 


2.77 


1.84 


1.08 


BM 


0.52 


0.88 


0.63 


0.64 


0.64 


0.64 


Anisotropic Huber-Ll 


0.59 


1.35 


0.57 


0.57 


0.57 


0.57 


BM + Anisotropic Huber-Ll 


0.25 


0.24 


0.26 


0.28 


0.21 


0.21 



Table A3. Position Error results after being evaluated with different methods. 

Methods Ml M2 M3 M4 M5 M6 

Lucas & Kanade 63.75 17.38 60.66 75.02 76.37 42.99 

BM+ Lucas & Kanade 53.35 45.51 57.63 39.68 29.86 19.55 

BM 1.30 1.49 1.80 1.81 1.49 1.50 

Anisotropic Huber-Ll 1.55 1.48 1.67 1.67 1.30 1.31 

BM + Anisotropic Huber-Ll 1.75 1.64 1.89 1.89 1.54 1.53 
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Table A4. Diameter error results after being evaluated with different methods. 
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M4 


M5 
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Lucas & Kanade 


4.36 


3.08 


3.79 


3.70 


2.52 
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BM+ Lucas & Kanade 


9.30 


5.79 


2.52 


2.25 


1.44 


1.02 


BM 


0.65 


0.86 


0.44 


0.48 


0.43 


0.44 


Anisotropic Huber-Ll 


0.47 


0.47 


0.37 


0.42 


0.38 


0.38 


BM + Anisotropic Huber-Ll 


0.26 


0.28 


0.26 
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0.20 


0.21 
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